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Abstract 

The Hubbard model is investigated in the framework of lattice density functional theory (LDFT). 
The single-particle density matrix jij with respect the lattice sites is considered as the basic 
variable of the many-body problem. A new approximation to the interaction-energy functional 
W[y] is proposed which is based on its scaling properties and which recovers exactly the limit of 
strong electron correlations at half-band filling. In this way, a more accurate description of W 
is obtained throughout the domain of representability of jij, including the crossover from weak 
to strong correlations. As examples of applications results are given for the ground-state energy, 
charge-excitation gap, and charge susceptibility of the Hubbard model in one-, two-, and three- 
dimensional lattices. The performance of the method is demonstrated by comparison with available 
exact solutions, with numerical calculations, and with LDFT using a simpler dimer ansatz for W. 
Goals and limitations of the different approximations are discussed. 

PACS numbers: 71.15.Mb, 71.10.Fd 
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I. INTRODUCTION 



Density functional theory (DFT) provides a rigorous framework for studying the physics 
of the many-body problem . 1 ! 2 ! 3 ' 4 The fundamental concept behind DFT is to replace the 
conventional dynamical variables that completely define a many-particle state (e.g., the 
wave function in a quantum mechanical problem, or the particle positions and momenta in 
a classical system) by considering the particle-density distribution p(f) as the basic variable. 
For this purpose the energy E of the system is expressed as a functional of p(r) separating 
the various energy terms in two main groups. The first one contains the contributions that 
depend explicitly on the problem under study, which result from the coupling between p(f) 
and the external fields V ext (r). The second one concerns the intrinsic energy of the many- 
particle system, namely, the kinetic energy T and the interaction energy W of the particles. 
These terms are universal functionals of p{r) in the sense that they are independent of the 
problem under study. While the general functional dependence of T[p] and W\p\ is not known 
explicitly, they can be formally expressed as the result of integrating out the microscopic 
degrees of freedom. In the case of ground-state electronic properties, this is achieved by 
imposing that T[p] and W[p] correspond to the minimum possible value of T + W for a 
given p(r) These basic notions have general validity and are therefore relevant to a wide 
variety of situations which may have very different physical origins. 

In the present paper the concepts of DFT are applied to investigate the physics of strongly 
correlated electrons in a narrow energy band. The theoretical description of these sys- 
tems is usually based on lattice Hamiltonians such as Anderson, 6 Hubbard,— Pariser-Parr- 
Pople, 8 and related models which focus on the most relevant electron dynamics at low 
energies . 9 ! 10 ! 11 ! 12 The study of many-body lattice models in the framework of DFT seems 
particularly interesting from various perspectives. On the one side, a detailed understand- 
ing of the electronic properties in the strongly correlated limit constitutes an important 
theoretical challenge. Exact results are rare or numerically very demanding and a variety 
of elaborate many-body techniques are specifically developed for their study. Therefore, a 
density functional approach with an appropriate ansatz for W should be a useful alternative 
tool for investigating at least some aspects of this complex problem. On the other side, one 
would like to extend the range of applicability of DFT to strongly correlated phenomena, like 
the separation of charge and spin degrees of freedom or the correlation induced localization, 
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where conventional local density approximations (LDAs) or generalized gradient approxi- 
mations (GGAs) are known to fail systematically. Moreover, the development of DFT on a 
lattice constitutes an intrinsically inhomogeneous approach, which provides a true alterna- 
tive to methods relying on the homogeneous electron gas. In this context, investigations on 
many-body models should open new insights into the properties of the interaction-energy 
functional that could also be useful for future extensions to more realistic Hamiltonians and 
first principles calculations. 

In past years, a number of density functional studies of lattice models have performed 
concerning in particular the determination of band-gaps in semiconductors,— the role of off- 
diagonal elements of the density matrix and the non-interacting v representability in strongly 
correlated systems,— or the development of energy functionals of the density matrix with 
applications to Hubbard and Anderson models.— In previous papers we have formulated a 
lattice density functional theory (LDFT) of many-body models by considering the density 
matrix 7^ with respect to the lattice sites % and j as the fundamental variable . 161 lls The 
interaction energy W of the Hubbard model has been calculated exactly as a function of 
jij for various periodic lattices having 7^ = 712 for nearest neighbors (NNs) i and j. On 
this basis, a simple general approximation to 1^(712) has been proposed which derives from 
exact dimer results, the scaling properties of W, and known limits. Using this ansatz, 
several ground-state properties of one-dimensional (ID) and two-dimensional (2D) systems 
have been obtained in good agreement with available exact solutions and accurate numerical 
calculations?^ In addition, applications to dimerized chains provided us with systematic 
results for the ground-state energy and charge excitation gap of the ID Hubbard model as a 
function of hopping alternation and Coulomb repulsion, including the crossover from weak 
to strong correlations.— LDFT appears therefore as an efficient method of determining the 
electronic properties of many-body lattice models, thus encouraging further developments 
and applications. 

The main purpose of this paper is to present a new approximation to the interaction- 
energy functional ^[7] of the Hubbard Hamiltonian and to apply it to determine several 
electronic properties of this model in the framework of LDFT. In Sec. |H]the basic formula- 
tion of LDFT is briefly reviewed. Different approximations to interaction-energy functional 
are presented and discussed in Sec. IIHI First, we analyze the properties of a previously 
proposed dimer ansatz^ and discuss its goals and limitations by comparison with known 
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exact results. Some shortcomings of this functional in the limit of strong correlations at 
half-band filling are pointed out. In order to overcome them, we propose a new approxima- 
tion based on the scaling properties of W, which recovers the correct behavior in the limit 
of strong interactions. In this way a more accurate description of the dependence on 7^ in 
different dimensions and lattice structures is obtained. The following sections are mainly 
concerned with applications to ID, 2D and 3D Hubbard models. Results for the ground- 
state energy, charge-excitation gap, and charge susceptibility are presented in Sees. IIVI IVl 
and IVH respectively. Comparison is made with the simpler dimer ansatz and with exact 
analytical or numerical solutions, whenever available, in order quantify the accuracy of the 
different approximations. Finally, Sec. I Vlll summarizes the conclusions and points out some 
perspectives of future developments. 

II. DENSITY-FUNCTIONAL THEORY OF LATTICE MODELS 

In order to be explicit we focus on the Hubbard model which is expected to capture the 
main physics of lattice fermions in a narrow energy band. The Hamiltonian 7 

H = tiAaCja + ftiifWr* (!) 

includes nearest neighbor (NN) hoppings t^, and on-site interactions given by U {h ia = 
clo-Cjo-). The hopping integrals ty are defined by the lattice structure and by the range of the 
single-particle hybridizations (typically, Uj = — t < for NN ij). They specify the system 
under study and thus play the role given in conventional DFT to the external potential 
Vext{r)- Consequently, the basic variable in LDFT is the single-particle density matrix 7^. 
The situation is similar to the density-matrix functional theory proposed by Gilbert for 
the study of nonlocal pseudopotentials,— since the hoppings are nonlocal in the sites. A 
formulation of DFT on a lattice in terms of the diagonal 7^ alone is possible only if one 
restricts oneself to models with constant for i ^ j. In this case the functional ^[7^] 
depends on the values of tij for i 7^ j and in particular on [//£.— 

The ground-state energy E gs and density-matrix 7^ are determined by minimizing the 
energy functional 

E[ 1 ]=E K [ 1 ] + W[ 1 ] (2) 
with respect to 7^. E[y] is defined for all density matrices that derive from a physical state, 
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i.e., that can be written as 



la = E^ = E<*l4c;«r|*> , (3) 

cr a 

where |\&) is an iV-particle state. Such 7^ are said to be pure-state iV-representable. An 
extension of the definition domain of E to ensemble-representable density matrices is 
straightforward following the work by Valone . 20 ^ 21 The first term in Eq. (J2J) is the kinetic 
energy associated with the electronic motion in the lattice. It is given by 

#*[7] =EVy,y, (4) 

ij 

thus including all single-particle contributions. The second term is the interaction-energy 
functional given by^ 

(r,) 



W[y) = min 

\Ef — >'y 



^E<*[7ll^l*[7]> 

i 

where the minimization runs over all iV-particles states ^[7]) that satisfy 

Efel%]) = Tii (6) 

a 

for all i and j. W[j] represents the minimum value of the interaction energy compatible 
with a given density matrix 7^. It is a universal functional of 7^ in the sense that it is 
independent of ty, i.e., of the system under study. However, W depends on the number 
of electrons N e , on the structure of the many-body Hilbert space, as given by N e and the 
number of orbitals or sites N a , and on the form of the model interactions.— 
E[y] is minimized by expressing 

7y = E 7ij<7 = E u ik*VkvU* jka (7) 

o" ka 

in terms of the eigenvalues r]k a (occupation numbers) and eigenvectors uq-u (natural orbitals) 
of Lagrange multipliers /1 and \k a (£fc<7 = ^ka/f]ka) are introduced in order to impose 
the constraints J2kaVka = N e and J2i \ u ika\ 2 = 1- Deriving with respect to u* ka and r\ ka 
(0 < T] ka < 1), one obtains the eigenvalue equations^^ 

/ , I Hj + — J Uika — ^ktjUjka , (o) 

with the subsidiary conditions £& CT < if f]k a = 1, £ka = /x if < rjka < 1, and Ska > /•* if 
rjka = 0. Self-consistency is implied by the dependence of dW/djij^ on r/ fccr and Ui\. a - This 



formulation is analogous to density-matrix functional theory in the continuum.— However, it 
differs from KS-like approaches which assume non-interacting f-representability and where 
only integer occupations are allowed . 13 i 14 In the present case, the fractional occupations of 
natural-orbitals play a central role. One may in fact show that in general < % CT < 1 for all 
ka. Exceptions are found in very special situations like the uncorrelated limit (U = 0) or the 
fully-polarized ferromagnetic state in the Hubbard model (S z = min{iV e , 2N a — N e }/2). This 
can be understood from perturbation-theory arguments — none of the rjkcr is a good quantum 
number for U ^ — and is explicitly verified by exact solutions of the Hubbard Hamiltonian 
on finite systems or the ID chain.— Therefore, all e^ a in Eq. (JHJ) must be degenerate and 
consequently the ground-state density matrix satisfies 

dW 

Uj + 77— = Oij n . (9) 

lij (J 

Notice the importance of the dependence of W on the off-diagonal density-matrix elements 
7ij which measure the degree of electron localization. Approximations of W in terms of 
the diagonal 7^ alone are not applicable in this framework (i^ 7^ for NN ij). Eq. © 
provides a self-consistent scheme to obtain the ground-state 7^ according to the variational 
principle. In the following section we present and discuss simple explicit approximations 
to W[j] that are intended to describe the electronic properties of the Hubbard model in 
different interaction regimes, band-fillings, and lattice structures. 



III. INTERACTION-ENERGY FUNCTIONAL IN THE HUBBARD MODEL 

The general functional W [7] , valid for all lattice structures and for all types of hybridiza- 
tions, can be simplified at the expense of universality if the hopping integrals are short 
ranged. For example, if only NN hoppings are considered, the kinetic energy Ek is in- 
dependent of the density-matrix elements between sites that are not NNs. Therefore, the 
constrained search in Eq. (jSJ) may be restricted to the ^[7]) that satisfy Eq. (jUJ) only for 
i — j and for NN ij. This reduces significantly the number of variables in W[j] and renders 
the determination and interpretation of the functional dependence far simpler. In particular 
for periodic lattices the ground-state jfj is translational invariant. Therefore, in order to 
determine E gs and jfj, one may set 7^ = n = N e /N a for all sites i, and 7^ = 712 for all NN 
pairs ij. In this case the interaction energy can be regarded as a simple function ^(712) 
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of the density-matrix element between NNs. This is certainly a great practical advantage. 
However, it should be noted that restricting the minimization constraints in Eqs. (0) and 
(JBJ) to NN also implies that W loses its universal character, since the NN map and the 
resulting dependence of W on 7 12 are in principle different for different lattice structures. 

The difficulties introduced by the lack of universality can be overcome by taking advantage 
of the scaling properties ^(712)- Recent numerical studies^ have in fact shown that W is 
nearly independent of system size N a , band filling n = N e /N a and lattice structure, if W is 
measured in units of the Hartree-Fock energy Z?hf — Un 2 /4 and if 712 is scaled within the 
relevant domain of representability [7^,712]- Here, 7° 2 stands for the largest possible value 
of the NN bond order 7 12 for a given N a , n, and lattice structure. It represents the maximum 
degree of electron derealization and corresponds to the uncorrelated limit. On the other 
side, 7^ refers to the strongly correlated limit of 712, i.e., to the largest NN bond order that 
can be obtained under the constraint of vanishing W. For half-band filling 7^ = 0, while 
for n 7^ 1, 7^ > 0.— Physically, the possibility of scaling the interaction energy means that 
the relative change in W associated to a given change in the degree of electron localization 
9\2 = (712 — 7i2)/(7i2 — 7i?) can be regarded as nearly independent of the system under 
study. This pseudo-universal behavior of W/Eft-p as a function of can be exploited to 
obtain good general approximations to 1^(712) by applying such a scaling to the functional 
dependence derived from a simple reference system or from known limits. 

In a previous paper we have proposed an approximation to the interaction energy W of 
the Hubbard model by extracting the functional dependence from the exact result for the 
Hubbard dimer, which is given byH 



This very simple expression satisfies several general properties of the exact 1^(712): 

(i) For 712 = 7i2, = E^f since the underlying electronic state ^[7i 2 ] is a single Slater 
determinant. Moreover, one observes that dW^ / dji2 = 00 for 7 12 = 7° 2 . This is a necessary 
condition in order that 7f 2 < 7° 2 already for arbitrary small U/t 7^ 0, as expected from 
perturbation theory. 

(ii) H /( - 2 - ) (7i 2 ) decreases monotonously with decreasing 712 reaching its lowest possible value, 
W = 0, for 712 = 7^. In other words, a reduction of the interaction energy is obtained at 
the expense of electron derealization. 




(10) 
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(iii) In the strongly correlated limit (7 12 <C 7° 2 ) one observes that oc ^\ 2 . Therefore, 
for U/t ^> 1, 7 9S oc t/U and E gs oc t 2 /U, a well known result in the Heisenberg limit of the 
half-filled Hubbard model.— 

A correct description of these basic properties and of the dependence of W/E HF on g 12 are 
at the origin of the remarkable performance of this simple dimer ansatz in the description 
of several ground-state properties of the Hubbard model.— 

In order to discuss the strongly correlated limit of Eq. (|1U|) in more detail we expand 
to lowest order in 712. At half-band filling one obtains 



with «2 = {li2)~ 2 - The exact interaction-energy W ex is also proportional to t/712 in the 
limit of small 712. Therefore, W ex can be expanded in the same form as Eq. (fTTj) but with 
a somewhat different coefficient a ex . Notice that in the case of the Hubbard dimer, we have 
7° 2 = 1 and af™ 1 = 1, which coincides of course with the exact result. Considering for 
example the ID chain, the 2D square lattice, and the 3D simple-cubic lattice one finds that 
the leading coefficients resulting from Eq. (fTT)|) are a\ = (vr/2) 2 ~ 2.47, a 2 = 6.09, and 
a-P = 9.30. These can be compared with the corresponding exact result derived from the 
Bethe-ansatz solution of the ID Hubbard chain,— or with perturbation-theory calculations 
for the square and simple-cubic lattices,— which are given by a l e ® = 2/ In 2 ~ 2.89, a 2 ® = 
6.91, and o?® = 10.94. One observes that Eq. (fTUj) reproduces correctly the trends in a with 
increasing dimensions. However, there is also a systematic underestimation of the interaction 
energy of the order of 12-15%. These quantitative discrepancies have direct consequences on 
the predicted properties, since the behavior of W for small 712 determines the ground-state 
density matrix 75^ and energy E gs in the strongly correlated limit. In fact, approximating 
W as in Eq. ()11|). writing the kinetic energy as = zt r ji 2l where z is the coordination 
number, and using Eq. (JHJ), one obtains 712 = (4:z/a)(t/U) and E gs = —{2z 2 /a){t 2 /U). 
Thus, an inaccuracy in a results in a similar relative error in 7^ and E gs for U/t^> 1. 

To overcome these shortcomings more flexible approximations to the interaction-energy 
functional are needed, which allow one to go beyond Eq. (|10p. Therefore, we propose a 
general ansatz of the form 



W<*> = (l/8)a 2 t/7i 2 2 + 0(^2) 



(11) 




(12) 
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where P n (gi2) is a function of gu = (712 — 712) /(7i2 — 7w)> thus incorporating the scaling 
properties of W without loss of generality. P n (<7i 2 ) is approximated by an n-order polynomial 
Pn(gn) = J2k=o a k9i2- This is justified by the fact that (W — E HF ) 2 is in general a well- 
behaved function of (712, even in the uncorrelated limit where dW/dj^ diverges (#12 = 1). 
The coefficients a& are to be determined from known properties of W . First of all, one 
observes that at half-band filling, and for bipartite lattices in general, the sign of 712 can be 
changed without altering W . Thus, P n {gv2) is an even function of #12 and = for odd k 
("fi2 = —^12 and 7^~ = — 7i2 + )- I n non-bipartite lattices away from half-band filling one 
may also set for simplicity = for odd k, since the dependence on is very similar for 
positive and negative 712, once the different domains of representability are scaled . 16 > 24 

The uncorrelated and fully-correlated limits of W (W = E HF for g 12 = 1, and W = 
for gu = 0) impose two simple conditions on the a^, namely, P n (l) = J2k a k = and 
P„(0) = a = 1. This defines the second-order approximation completely. In this case, 
^2(^12) = 1 — #12; which coincides with the above discussed dimer ansatz [Eq. (fT0|) ]. The 
two approaches are therefore consistent. The dimer ansatz can be regarded as the simplest 
polynomial-based approximation of W, as given by Eq. (|12jl . that satisfies the obvious limits. 

The 4th-order approximation introduces the aimed additional flexibility that can 
be exploited to reproduce the strongly correlated limit of W exactly. Expanding Eq. (JT2J) 
to second-order in 7 12 one observes that at half-band filling this is achieved when a = 1, 
a>2 = — a e x(7i 2 ) 2 = — a ex/&2, and 0,4 = — (ao + a 2 ). Thus, the 4th-order approximation to W 



where k = a ex / 02 > is the ratio between the small-712 expansion coefficients of W ex and 
W( 2 \ The value of n depends on the lattice structure or system dimensions. At half-band 
filling it can be determined by applying perturbation-theory to the Heisenberg limit of the 
Hubbard model. 25 For instance, for the ID chain, 2D square lattice, and 3D simple-cubic 
lattice one obtains, respectively, k\d = 8/(7r 2 ln2) = 1.169, k 2 d = 1.135, and k S £, = 1.176. 
Notice that k depends rather weakly on the lattice structure and that it is not very far from 
the dimer value k = 1, for which Eq. ()13j) reduces to Eq. (fTTH) . Therefore, the 4th-order 
term appears as a relatively small correction to the 2nd-order approximation. Higher-order 
polynomial approximations to W could be derived in an analogous way, provided that reliable 
information is available on the following terms of the small-7jj expansion of W ex . This gives 



is given by 




(13) 
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Y-12^12 

FIG. 1: Interaction energy W(ji2) of the one-dimensional (ID) Hubbard model at half-band filling 
(n = 1) as a function of the density-matrix element or bond order 712 between nearest neighbors. 
7i 2 refers to the ground-state bond order in the uncorrelated limit (U = 0). Results are given 
for the dimer approximation [Eq. (|lUj) . dashed], the 4th-order approximation [Eq. P|l 

with k = k\d = 1.169, solid], and the exact W ex [Eq. Q, crosses] which is derived from the Bethe 
ansatz solution^ In the inset the corresponding relative errors are shown. 

the possibility of further improving the accuracy of the results by incorporating a more 
detailed description of the strongly correlated limit. 

In Fig. n the approximate interaction energies and of the half-filled ID Hub- 
bard chain are compared with the corresponding exact result W ex , as derived from the 
Bethe- Ansatz solution.— As already observed,— even the simplest dimer ansatz follows 
W / ex(7i2) quite closely all along the crossover from weak to strong correlations. In this case 
the interaction energy is always underestimated, and the absolute value of the relative er- 
ror e = \W — W ex \/W ex increases monotonously as 712 decreases, reaching about 15% for 
712/712 < 0.4. The 4th-order approximation provides a significant advance, not only for 
712/712 <C 1 but in the complete domain of representability. For the relative error e is 
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reduced to less than 1% for 712/712 < 0.4 (e — > for 712 — > 0). The largest discrepancies 
are found for 712/712 — 0.8-0.9, where e reaches only 3%. An appreciable improvement in 
the accuracy of the derived properties can be therefore expected. In the following sections, 
Eqs. (jlUj) and (J13j) are applied in the framework of LDFT to determine several electronic 
properties of the Hubbard model in ID, 2D, and 3D periodic lattices. 

IV. GROUND STATE ENERGY 

In Fig. |21 the ground-state energy E gs of the half-filled ID Hubbard model is given as a 
function of the Coulomb repulsion strength U/t. Comparison between LDFT and the Bethe- 
Ansatz exact solution shows that 4th-order approximation improves significantly the already 
good results derived using the dimer ansatz.— This concerns not only the strongly correlated 
limit where, as expected, recovers the exact result, but the complete range of U/t. The 
largest quantitative discrepancies between exact and 4th-order results are in fact very small. 
They amount to less than 3% and are found for intermediate interaction strengths (U/t ~ 4). 
In contrast, the relative error in the dimer ansatz increases monotonously with U/t reaching 
about 17% for U/t = 00 (see the inset of Fig. It is interesting to note that in both 
cases no artificial symmetry breaking is required in order to describe correctly the electron 
localization induced by correlations and the resulting dependence of E gs on U/t, as it is 
often the case in other mean-field approaches. 

The higher performance obtained with the 4th-order correction originates in an improved 
accuracy of both kinetic and Coulomb contributions to the ground-state energy. As shown 
in Fig. El the kinetic energy Ek < increases monotonously with increasing U/t, first rather 
slowly up to U/t c± 4, and then more rapidly when electron localization starts to set in. For 
U/t < 4, the values of E K obtained using W^ 2 > and W^> are very close to the exact result 
(typically \E ( p - Ef\/Ef < 2.6% and \E$ - Ef\/Ef < 2.0%). For U/t > 4 the dimer 
ansatz shows some limitations while the 4th-order approximation remains very accurate (for 
example, for U/t = 12, \eP - Ef\/ Ef ~ 13% and \E$-E%\/E% < 2.4%). The Coulomb 
energy Ec shows the usual non- monotonous behavior, first increasing with U/t in the weakly 
correlated regime and then decreasing as the strongly-correlated limit is approached. This 
behavior is correctly described by both 2nd- and 4th-order approximations. However, one 
finds that it is in general more difficult to accurately describe Ec as compared to Ek- The 



11 




0.0 0.2 0.4 0.6 0.8 1.0 



U/(U+4t) 

FIG. 2: Ground-state energy E gs of the half-filled ID Hubbard model as a function of the Coulomb 
repulsion strength U/t. The dashed curves refer to lattice density-functional theory (LDFT) using 
the dimer approximation to W [Eq. (|10|) ] and the solid curves to the 4th-order approximation 
[Eq. with k = k\d = 1.169]. The crosses are the exact results derived from the Bethe-ansatz 
solution^ The corresponding relative errors are given in the inset. 

2nd-order E c underestimates (overestimates) E(f appreciably for 2 < U/t < 5 (U/t > 5). 
The 4th-order correction provides a clear improvement over the dimer ansatz, by increasing 
Ec in one case (U/t < 5) and reducing it in the other (U/t > 10). As for E gs , the remain- 
ing differences with the exact results are quite small and correspond to intermediate U/t. 
Summarizing, one may observe that the accuracy of the calculated E gs is not the result of 
a strong compensation of errors, since a very good performance is achieved for the kinetic 
and Coulomb energies separately. 

In Fig. |3] results are given for E gs of the 2D square lattice and 3D simple cubic lattice 
at half-band filling. For U/t < 3 the 2nd-order and 4th-order results are almost indistin- 
guishable, while for U/t > 4 the 4th-order approximation yields somewhat higher values 
(Ef) < E$ < 0). These trends are very similar to those observed in the ID chain. The 
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U/(U+4t) 

FIG. 3: (a) Kinetic energy Ek and (b) Coulomb energy Ec of the half-filled ID Hubbard model 
as a function of U/t. The dashed curves correspond to the dimer approximation [Eq. (|1U[)]. the 
solid curves to the 4th-order approximation [Eq. (|13j) with k = k\d = 1.169], and the crosses to 
the Bethe-ansatz exact solution^ 

LDFT calculations for 2D and 3D systems compare well with far more demanding quantum 
Monte Carlo (QMC) studies^*^ (see Fig. EJ). Furthermore, the reliability of the LDFT re- 
sults is confirmed by comparison with exact Lanczos diagonalizations on small clusters, for 
example, on a N a = 3 x 4 cluster of the square lattice with periodic boundary conditions. 
In this case, like in ID, the overall performance is very good, with the largest quantitative 
discrepancies being observed for intermediate values of U/t. For instance, for U/t = 4 one 
obtains \E$-E%\/\E%\ = 4.2 x KT 2 , and for U/t = 16 \E$-E%\/\E%\ = 3.2 x 1CT 2 . In 
conclusion, LDFT using Eq. (fTSj) for the interaction energy W yields an accurate description 
of the ground-state energy of the Hubbard model in different dimensions.— 
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FIG. 4: Ground-state energy E gs of the half-filled Hubbard model as a function of the Coulomb 
repulsion U/t: (a) two-dimensional (2D) square lattice and (b) three-dimensional (3D) simple cubic 
lattice. The dashed curves correspond to the dimer approximation [Eq. (jlOj) ] and the solid curves 
to the 4th-order approximation [Eq. ()13|) with (a) k = k>2D = 1-135 and (b) k = k^d = 1-176]. The 
crosses with error bars refer to quantum Monte Carlo (QMC) calculations 

V. CHARGE EXCITATION GAP 

The charge-excitation or band gap 



is a property of considerable interest, which measures the low-energy excitations associated 
to changes in the number electrons N e , and which is very sensitive to the degree of electronic 
correlations. It can be related to the discontinuity in the derivative of the ground-state 
kinetic energy E K and correlation energy E corT = E c — -EW with respect to band-filling n. 
For N a — > oo and n = 1, it is given by 




(14) 




(15) 
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FIG. 5: Charge excitation gap AE C of the Hubbard model at half-band filling as a function of U/t: 
(a) ID chain, (b) 2D square lattice, and (c) 3D simple cubic lattice. The dashed curves correspond 
to the dimer approximation [Eq. (|lUj) ] and the solid curves to the 4th-order approximation [Eq. Q13|)] 
with (a) k = kid = 1-169, (b) k = nop = 1.135, and (c) n = k^d = 1-176. The crosses refer to 
exact results in the ID chain (Ref. 23J) and to QMC calculations in 2D (Ref. \2(i \ and 3D lattices 



(Ref. 



23). 



where e = [Ek + E corr )/N a . The determination of AE C is in general a more difficult task 
than the calculation of ground-state properties like E gs , Ek, and Eq. In fact, the band-gap 
in semiconductors has been an important problem which motivated numerous works in the 
context of DFT in the continuum. Therefore, AE C appears as a particularly interesting 
property to investigate with the present lattice density-functional formalism. 

In the half-filled Hubbard model on bipartite lattices, AE C increases with increasing U/t 
(AE C — for U/t — 0) and approaches the limit AE C —>([/ — Wb) for U/t — ► oo, where 
Wb is the width of the single-particle band (wb = 4t, 8t, and 12t for the ID, 2D square, 
and 3D simple-cubic lattices, respectively). Fig. El presents LDFT results for AE C in ID, 
2D, and 3D Hubbard models (n = 1). Comparison with the exact Bethe-Ansatz solution 
for the ID chain^ and with available QMC calculations for the square^ and simple cubio^l 
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lattices shows a good overall agreement. It should be however noted that in the ID case 
the gap is significantly overestimated for U/t <C 1. Here we obtain AE C oc (U/t) 2 , while in 
the exact solution AE C increases much more slowly, namely, exponentially in —t/U. This 
discrepancy concerns both the 2nd-order and the 4th-order approximations, which are nearly 
indistinguishable for U/t < 2-4. Consequently, it is possible that the results for 2D and 3D 
lattices shown in Fig. |3] also overestimate the gap for small U / 1. In any case, the accuracy 
of LDFT improves rapidly with increasing U / 1, as electron localization starts to set in, and 
the error in AE C tends to zero for large U/t. Therefore, the development of a Mott insulator 
in the strongly correlated limit is correctly described. 

It is important to remark, in the context of metal-insulator transitions in three dimen- 
sions, that our calculations on the SC lattice yield a finite gap AE C > for n — 1 and 
all U/t > 0. This is consistent with previous results on 3D bipartite lattices, which are 
expected to be antiferromagnetic (AF) insulators for all U/t > . n i 12 The functionals 
and reproduce correctly this behavior, as well as the formation of local moments 

(Sf) = 3(1 — 2(nj|nj|))/4, without involving a spin-density-wave symmetry breaking. This 
can be understood be recalling that they are based on the exact functional of the Hubbard 
dimer which, being a bipartite cluster, incorporates AF correlations (n = 1). However, the 
properties change qualitatively if frustrations become important (e.g., in non-bipartite lat- 
tices or if second NN hopping are significantly large). In this case it has been shown that 
the half-filled Hubbard model is a metal with AE C = for small U > and that a metal- 
insulator transition takes place at a finite interaction strength U c , which is of the order of 
the single-particle band width Wb-— This behavior is not reproduced by the functionals 
and W {i \ even if they are applied to compact lattices (e.g., the face-centered cubic lattice), 
since they are free from any singularities throughout the domain of representability (except 
for gi2 = 1) and since the resulting 7?? are smooth functions of U/t. Notice that the exact 
functional W ex may show a far more complex behavior, particularly if the nature of the 
state I ^[7]) yielding the minimum of Eq. © changes as a function of 7. This is expected 
to be the case at a metal- insulator transition, where a discontinuous decrease of {h^fiii) 
occurs. Finally, let us point out that the results derived from Eqs. (fHH) and ()13|) for large 
U/t (U > U c ) are consistent with previous studies. This concerns not only the presence of 
a finite gap AE C , but also the fact that the number of double occupations does not vanish 
on the insulating side of the transition . n ' 12 
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Comparing 2nd- and 4th-order approximations for U/t > 2-4 one observes that the charge 
gap is always somewhat smaller in the latter case. For the ID chain, the reduction of AE C 
due to the 4th-order correction improves the agreement with the exact solution appreciably 
(e.g., |A£W - AE™\ ~ \AE^ - AE™\/2 for U/t > 10). In the considered 2D and 3D 
lattices, the differences between 2nd- and 4th-order results are similar to those observed in 
the ID chain. Comparison with QMC calculations shows a good overall agreement although 
some quantitative differences can be noted. For example, as shown in Fig. our values 
for AE C are somewhat smaller than the QMC ones for the 2D (3D) lattice with U/t — 4 
(U/t = 8) and somewhat larger for U/t = 8 {U/t = 12). In summary, the ensemble of ID, 2D 
and 3D results shows that LDFT provides a very simple and efficient method of calculating 
the charge excitation energies of the Hubbard model in different dimensions and interaction 
regimes. However, the proposed approximations to W are still not quite satisfactory in the 
weakly-correlated limit and deserve to be improved. 



VI. CHARGE SUSCEPTIBILITY 



The charge susceptibility Xc is defined by 

< 16 > 

where n = N e /N a is the number of electrons per site and fi the chemical potential. It 
represents the many-body density of electronic states at the Fermi energy \i and thus provides 
very useful information on the low-energy charge-excitation spectrum as a function of band 
filling. In Figs. EHE1 Xc is given as a function of fi for ID, 2D, and 3D Hubbard models 
on bipartite lattices for representative values of U/t. The LDFT calculations reported in 
these figures were performed using the dimer ansatz for W given by Eq. (jlOJ) . As it will 
be discussed below, the 4th-order approximation [Eq. ()13|)] yields very similar results. In 
the case of the ID chain comparison is made with the exact Xcin), which is obtained from 
the Bethe ansatz solution.^ For the 2D square lattice, we also show ground-state QMC 
results^ for U/t = 4. Notice that in bipartite lattices, as those considered here, electron- 
hole symmetry implies that \ c is the same for band fillings n and n! = 2 — n, and therefore 

Xc{v) = Xc(v' = U-fi). 

In the absence of interactions \ c coincides with the single-particle density of states of the 
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FIG. 6: Charge susceptibility Xc of the ID Hubbard model as a function of the chemical potential 
fx for different Coulomb repulsions U/t. The solid curves refer to LDFT in the dimer approximation 
[Eq. (jlUj) ] and the dashed curves with crosses to the exact Bethe ansatz solution^ For U/t = 64 
only the lower Hubbard band is shown. 

corresponding lattices (U = 0). These are gapless and show the usual van Hove singularities 
at the band edges fi = ^Wb/Z and at some points within the bands of the square and simple- 
cubic lattices. For finite U a gap AE C = n(n = 1 + ) — fi(n = 1~) opens at half-band filling 
which increases monotonously with U, as discussed in the previous section. Thus, the two 
so-called lower and upper Hubbard bands start to be distinguished, which correspond to 
hole and electron doping respectively. The separation of the bands becomes particularly 
clear for U ~ w^, when AE C reaches values of the order of single-particle band width Wb (see 
Figs.lBHSl)- At the same time the width of the lower and upper bands increases with U, from 
Wb/2 for U = + , to Wb for U = +oo. These qualitative features are common to bipartite 
lattices in all dimensions and are correctly described by LDFT. 

In the ID case, where a detailed comparison with the exact solution is possible, we observe 
that our results are very accurate except close to half-band filling and for small or moderate 
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FIG. 7: Charge susceptibility Xc of the Hubbard model on a 2D square lattice as a function of 
the chemical potential /x for different Coulomb repulsions U/t. The solid curves are obtained using 
LDFT and the dimer approximation to the interaction-energy functional [Eq. Q1U|)]. The crosses for 
U/t = 4 refer to ground-state QMC calculations (Ref. l26h. For U/t = 64 only the lower Hubbard 
band is shown. 

values of U/t (see Fig. El). The nature of the discrepancies close to n = 1 is basically twofold. 
First, we find again the overestimation of the band gap AE C which is relatively important 
for small U/t (see also Sec. fV)l . Consequently, the band edges //(n = 1~) and /z(n = 1 + ) 
are not precisely reproduced in this limit, even if the absolute error e = |/i — fi ex \ always 
remains reasonably small. The largest inaccuracies are found for U/t ~ 3 and amount to 
e/wb = 8.1 x 10~ 2 . Nevertheless, this problem disappears as U/t increases, since e tends 
rapidly to zero in the strongly correlated limit (e.g., e/wb = 1.2 x 10 -4 for U/t = 16). The 
second limitation concerns the shape of \ c close to half-band filling. The exact solution of 
the ID chain shows sharp divergences in \ c at the gap edges fi(n = 1~) and fi(n = 1 + ) for 
U > 0, which we fail to reproduce. For small and moderate U/t, for example U/t = 1 or 
4, we obtain a nearly constant Xc for A* ~^ M ra = l^) , while the exact result is Xc +o° 
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FIG. 8: Charge susceptibility Xc of the Hubbard model on a 3D simple cubic lattice as a function 
of the chemical potential [i for different Coulomb repulsions U/t. The results are obtained using 
the dimer approximation to the interaction-energy functional [Eq. (jlUj) ]. For U/t = 64 only the 
lower Hubbard band is shown. 

(see Fig. EJ). Notice, however, that the increase and divergence of xT are sharply localized 
in a narrow range of particularly for small U/t. The divergence of Xc for n — > 1 could be 
reproduced by considering broken symmetry solutions of the LDFT equations, like in the AF 
Hartree-Fock approximation. Even so, it would be more interesting to describe this effect 
without involving a symmetry breaking, which is known to be artificial, and which could 
affect the results on the kinetic, Coulomb, and total energies, particularly in the case of finite 
systems.— As we approach the strongly correlated limit the LDFT results for Xc develop 
peaks at the gap edges, which height increases with U/t, thus approaching asymptotically 
the exact result. Still, Xc always remains finite for all finite U/t (see Fig. El for U/t = 16 
and 64). The very good performance for large U/t can be understood by recalling that for 
U/t = +oo, the LDFT results correspond to the fully-polarized or Nagaoka state,— which is 
the exact ground-state in ID for all n (U/t = +oo).— In higher dimensions it is possible that 
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our calculations yield a finite-height peak for n — > 1, where a true divergence of \ c could be 
present. This seems to be the case in the 2D square lattice where we observe narrow peaks 
in Xc at the band edges. In fact ground-state QMC calculations on the square lattice with 
U/t = 4 predict a divergent \ c at half-band filling (see Fig. |7j). In contrast, no such peaks 
are found in our calculations of Xc for the 3D simple-cubic lattice (see Fig. [SJ). 

As already discussed in the previous section, it is important to remark that the results 
presented in Fig. |H] are representative of bipartite lattices which at the half-band filling 
show an AF insulating behavior for all U/t > 0. In this case of our results are in good 
qualitative agreement with previous studies. The obtained simple Hubbard-approximation- 
like structure of Xo with a lower and upper Hubbard bands, also applies to frustrated 
lattices or to paramagnetic phases provided that U/t is sufficiently large to bring the system 
on the insulating side of the metal-insulator transition (U > Ue) . 11 ^ 12 However, it has been 
shown that the presence of frustrations drives the system into an AF metallic state at small 
U/t which contrasts with the AF insulator found in the absence of frustrations.— In this 
case the spectral density presents — in addition to a progressive development of lower and 
upper Hubbard bands with increasing U/t — a Kondo-like resonance at half-band filling 
(/i — U/2 = 0) characterized by a constant-height peak having a width that decreases with 
increasing U/t and that vanishes at U c , i.e., at the transition to the insulating state (see 



Ref. 



111 ). The functionals or fail to reproduce this kind of behavior, even when 



applied to compact structures (e.g., fee lattice). This limitation does not seem surprising, 
since Eqs. (1TUJ) and (fT3|) were derived from the properties of a bipartite system, and since the 
extensions presented in this paper, while achieving an accurate W in the strongly correlated 
limit at half-band filling, do not aim a precise description at small U/t and as a function 
of n. Exploring the functional dependence of W in frustrated structures, particularly close 
to n = 1 and g±2 — 1, should provide very useful clues in view of developing practical 
approximations capable of describing these remarkable effects. 

Let us finally point out that we have also determined Xc using as approximate 

interaction energy [Eq. (fT3*|) ] and found that the results are very similar to those obtained 
using and shown in Figs. EH3 in both cases the results are extremely good away 

from half-band filling, nearly indistinguishable from the ID exact solution. Close to n = 1, 
the 4th-order calculations yield smaller AE C and thus perform slightly better for U/t < 4. 
However, the divergences of Xc at the gap edges are not reproduced. Therefore, the 4th- 
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order corrections do not provide a significant improvement over the dimer ansatz concerning 
Xc of the ID chain. This is probably related to the simple form considered for which 
uses a coefficient k that is independent of n [see Eq. (113)) ]. While this approximation seems 
satisfactory for applications that concern a fixed band filling, it appears as a limitation for 
properties like AE C or Xc, where a precise description of the dependence of W on n is crucial. 
For 2D and 3D lattices the 4th-order results for Xc are also very similar to those shown in 
Figs. Eland El 

VII. DISCUSSION 

A new approximation to the interaction-energy functional W[y] of the Hubbard model 
has been proposed in the framework of lattice density functional theory, which exactly 
recovers the limit of strong electron correlations at half-band filling. The simpler ansatz 
which was derived from the functional dependence of W in the Hubbard dimei^ is thereby 
extended and improved. A more accurate description of W is achieved throughout the 
domain of representability of 7^ including the crossover from weak to strong correlations. 
Several properties have been determined by applying this functional to one-, two-, and three- 
dimensional lattices. Ground state energies, as well as kinetic and Coulomb energies, were 
successfully determined in all dimensions and interaction regimes. Very good results are 
also obtained concerning the charge-excitation gap AE C and the charge susceptibility Xc of 
bipartite lattices, except very close to half-band filling (n = 1) and for small values of the 
Coulomb repulsion strength (U/t < 4). This reveals some limitations in the description of 
the band-filling dependence of W for n ~ 1 and 712 — 7i 2 , which deserve more detailed 
investigations. Further insight on the origin of this problem could be obtained, for example, 
by analyzing the properties of the exact W as derived from the Bethe ansatz exact solution of 
the ID chain and from Lanczos diagonalizations in finite 2D clusters with periodic boundary 
conditions. Moreover, the functional dependence of W could be determined in the limit of 
large 7 12 (i.e., 712 — > 7° 2 corresponding to the weak correlations) by applying perturbation 
theory for small U/t. In this way, more accurate approximations to W could be developed 
in order to improve the results on AE C and Xc in this limit, particularly concerning the 
differences between bipartite and non-bipartite lattices. 

Besides these methodological aspects, the accuracy of the results and the simplicity of 
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the calculations encourage new applications of the present approach to related problems of 
current interest like dimerized one- dimensional chains and ladders, the 2D square lattice with 
competing first and second nearest-neighbor hoppings, the properties of n electrons in doped 
fullerenes and nanotubes in the framework of Hubbard or PPP models, or the connection 
with the continuum's DFT using minimal basis sets. In this way, a novel density-functional 
route to the physics of strongly correlated fermions is opened. 
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